9. Γεωεπεξεργασία διανυσματικών δεδομένων¶

Ο χώρος και οι γεωμετρικές δομές του μπορούν να αναπαρασταθούν μέσω διανυσματικών δεδομένων (Vector). Οι τρείς βασικές δομές δεδομένων είναι:

Τα σημεία
Οι γραμμές
Τα πολύγωνα

alt text

Επιπλέον αυτά τα γεωμετρικά δεδομένα συνοδεύονται και απο περιγραφικά δεδομένα που αφορούν τις ιδιότητες ή τα χαρακτηριστικά αυτών των δεδομένων.

alt text

Στα Σ.Γ.Π. ο πιο συνηθισμένος τύπος αρχείων αποθήκευσης αυτών των δεδομένων είναι το shapefile. Πλέον έχουν αναπτυχθεί και άλλοι τύποι όπως geojson, geopackage και χωρικές βάσεις (geodatabases) όπως η Postgresql/Postgis.

Ο προσδιορισμός της γεωγραφικής θέσης στην υδρόγειο γίνεται μέσω ενός ζεύγους γεωγραφικών συντεταγμένων. Κάθε σημείο στον χώρο προσδιορίζεται γεωγραφικά από το γεωγραφικό μήκος (λ) και το γεωγραφικό πλάτος (φ).

alt text

Για να αποδοθεί η τρισδιάστατη υδρόγειος σφαίρα σε ένα δυσδιάστατο σύστημα αναφορά χρησιμοποιείται ένα προβολικό σύστημα.

alt text

Κάθε προβολικό σύστημα της σφαίρας στο επίπεδο εισάγει μια σειρά παραμορφώσεων που αφορά το σχήμα των γεωμετρικών δομών, την κλίμακα, την έκταση και τις αποστάσεις. Ανάλογα το προβολικό σύστημα κάποιες από τις παραπάνω παραμορφώσεις εμφανίζονται σε μεγάλο βαθμό και άλλες όχι. Οπότε ανάλογα το είδος της έρευνας ο ερευνητής οφείλει να γνωρίζειπαραμορφώσεις εισάγει η κάθε προβολή και ανάλογα να επιλέγει την προβολή με τα λιγότερα σφάλματα.

Python βιβλιοθήκες για διανυσματικά δεδομένα¶

Για την ανάγνωση, εγγραφή, επεξεργασία διανυσματικών δεδομένων στην Python έχουν καθιερωθεί μια σειρά βιβλιοθηκών. Η αρχαιότερη και βασική βιβλιοθήκη είναι η GDAL/OGR. Επειδή η βιβλιοθήκη δεν είναι ιδιαίτερα συμβατή σε σχέση με τον τρόπο συγγραφής της Python και επιρρεπής σε σφάλματα έχουν αναπτυχθεί πιο σύγχρονες βιβλιοθήκες όπως η βιβλιοθήκη [Fiona] που είναι ιδιαίτερα χρήσιμη για την ανάγνωση/εγγραφή διανυσματικών δεδομένων Fiona και η βιβλιοθήκη Shapely η οποία χρησιμοποιείται για την επεξεργασία και ανάλυση.

Η βιβλιοθήκη Shapely¶

shapely from WKT¶

Μπορούμε να δημιουργήσουμε αντικείμενα shapely που αναπαριστούν σημεία ή γραμμές ή πολύγωνα μέσω WKT. Η γλώσσα (WKT)[https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry] είναι μια ειδική διάλεκτος για την περιγραφή διανυσματικών αντικειμένων. Εισάγουμε τις απαραίτητες υπο-βιβλιοθήκες (geometry και wkt) από την βιβλιοθήκη shapely.

import shapely.geometry
import shapely.wkt

Πολύγωνα¶

Καλούμε την μέθοδο shapely.wkt.loads() για να δημιουργήσουμε shapely objects από wkt

pol1 = shapely.wkt.loads("POLYGON ((0 0, 0 -1, 7.5 -1, 7.5 0,0 0))")

pol1

Η εκτύπωση του shapely αντικειμένου επιστρέφει μία περιγραφή σε μορφή WKT.

print(pol1)

POLYGON ((0 0, 0 -1, 7.5 -1, 7.5 0, 0 0))

Εκτύπωση του τύπου του αντικειμένου pol1

type(pol1)

shapely.geometry.polygon.Polygon

pol2 = shapely.wkt.loads("POLYGON ((0 1, 1 0, 2 0.5, 3 0, 4 0, 5 0.5, 6 -0.5, 7 -0.5, 7 1, 0 1))")
print(pol2)

POLYGON ((0 1, 1 0, 2 0.5, 3 0, 4 0, 5 0.5, 6 -0.5, 7 -0.5, 7 1, 0 1))

pol2

Δημιουργία MultiPolygon shapely object από αντίστοιχη συλλογή πολυγώνων MULTIPOLYGON

pol3 = shapely.wkt.loads("""
MULTIPOLYGON
(((40 40, 20 45, 45 30, 40 40)),
((20 35, 10 30, 10 10, 30 5, 45 20, 20 35), (30 20, 20 15, 20 25, 30 20)))
""")

type(pol3)

shapely.geometry.multipolygon.MultiPolygon

pol3

type(pol3)

shapely.geometry.multipolygon.MultiPolygon

Σημεία¶

pnt1 = shapely.wkt.loads("""
POINT (30 10)
""")

pnt1

Συλλογή σημείων

pnt2 = shapely.wkt.loads("""
MULTIPOINT(0 0,1 1)
""")

pnt2

Γραμμές¶

line1 =shapely.wkt.loads("""
LINESTRING(1.5 2.45,3.21 4)
""")

line1

Συλλογή γραμμών

line2 =shapely.wkt.loads("""
MULTILINESTRING((0 0,-1 -2,-3 -4, -3 -8),(2 3,3 4,6 7))
""")

line2

Shapely μέσω συναρτήσεων¶

from shapely.geometry import Point, LineString, Polygon, MultiPoint, MultiLineString, MultiPolygon

Κάθε συνάρτηση λαμβάνει διαφορετικά ορίσματα ανάλογα με την συνάρτηση.

Στην συνέχεια ακολουθεί η δημιουργία ενός αντικειμένου shapely γεωμετρία σημείου.

pnt1 = shapely.geometry.Point((2, 0.5))
pnt1

Για την δημιουργία συλλογής σημείων multipoint δίνουμε σαν όρισμα στην συνάρτηση shapely.geometry.MultiPoint μία λίστα από πλειάδες (tuples). Η κάθε πλειάδα (x,y)αντιστοιχεί στις συντεταγμένες ενός σημείου.

coords = [(2, 0.5), (1, 1), (-1, 0), (1, 0)]
pnt2 = shapely.geometry.MultiPoint(coords)
pnt2

Για την δημιουργία γραμμών χρησιμοποιείται πάλι μία λίστα από tuples που περιγράφουν τις κορυφές της γραμμής. Στο παρακάτω παράδειγμα χρησιμοποιούμαι την προηγούμενη πλειάδα αλλά πλέον χρησιμοποιούμε την μέθοδο shapely.geometry.LineString για την δημιουργία γραμμής και όχι την μέθοδο hapely.geometry.MultiPoint που δημιουργεί συλλογές σημείων.

line1 = shapely.geometry.LineString(coords)
line1

Κατασκεύη γραμμής από Point Objects:

p1 = shapely.wkt.loads("POINT (0 0)")
p2 = shapely.wkt.loads("POINT (1 1)")
p3 = shapely.wkt.loads("POINT (2 -1)")
p4 = shapely.wkt.loads("POINT (2.5 2)")
p5 = shapely.wkt.loads("POINT (1 -1)")

shapely.geometry.LineString([p1, p2, p3, p4, p5])

Αντίστοιχα μπορούμε να δημιουργήσουμε συλλογές γραμμών.

Δημιουργούμε ξεχωριστά αντικείμενα γραμμών shapely και στην συνέχεια καλούμε την συνάρτηση shapely.geometry.MultiLineString όπου ορίζουμε σαν όρισμα μια λίστα με τις μεμονωμένες γραμμές.

l1 = shapely.geometry.LineString([(2, 0.5), (1, 1), (-1, 0), (1, -1)])
l2 = shapely.geometry.LineString([(-2, 1), (2, -0.5), (3, 1)])
line2 = shapely.geometry.MultiLineString([l1, l2])
line2

Αντίστοιχα δημιουργούμε πολυγωνικά αντικείμενα μέσω μια λίστας με πλειάδες σημείων που χρησιμοποιείται σαν όρισμα στην συνάρτηση shapely.geometry.Polygon

coords = [(0, 0), (0, -1), (7.5, -1), [7.5, 0], (0, 0)]
shapely.geometry.Polygon(coords)

Αν θέλουμε μπορούμε να περάσουμε μία δεύτερη λίστα η οποία περιλαμβανει επιμέρους λίστες με πλειάδες σημείων τα οποία περιγράφουν τρύπες μέσα στο πολύγωνο.

coords_exterior = [(0, 0), (0, -1), (7.5, -1), (7.5, 0), (0, 0)]
coords_interiors = [[(0.4, -0.3), (5, -0.3), (5, -0.7), (0.4, -0.7), (0.4, -0.7)]]
shapely.geometry.Polygon(coords_exterior, coords_interiors)

Με την συνάρτηση shapely.geometry.MultiPolygon δημιουργούμε αντίστοιχα συλλογές πολυγώνων από μεμονωμένα αντικείμενα shapely.geometry.polygon.Polygon

multipolygon1 = shapely.geometry.MultiPolygon([pol2, pol3])
multipolygon1

/home/leonidas/anaconda3/envs/book/lib/python3.10/site-packages/shapely/geometry/multipolygon.py:202: ShapelyDeprecationWarning: __getitem__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  shell = ob[0]
/home/leonidas/anaconda3/envs/book/lib/python3.10/site-packages/shapely/geometry/multipolygon.py:203: ShapelyDeprecationWarning: __getitem__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  holes = ob[1]

print(multipolygon1)

MULTIPOLYGON (((0 1, 1 0, 2 0.5, 3 0, 4 0, 5 0.5, 6 -0.5, 7 -0.5, 7 1, 0 1)), ((40 40, 20 45, 45 30, 40 40)))

Σύμφωνα με τις προδιαγραφές simple features το παραπάνω object δεν είναι έγκυρο γιατί ένα πολύγωνο τέμνει ένα άλλο σε άπειρο αριθμό σημείων. Μπορούμε να ελέγξουμε την εγκυρότητα ενός αντικειμένου με την κλήση τις ιδιότητας is_valid. Επίσης κατά την οπτικοποίηση στο προηγούμενο βήμα του αντικειμένου multipolygon1 αυτό εμφανίζεται με κόκκινο (αντί για πράσινο). Ένδειξη ότι δεν ειναι valid.

multipolygon1.is_valid

True

Ταυτόχρονα μπορούμε να φτιάξουμε σύνθετες συλλογές από επιμέρους αντικείμενα shapely.

geo_collection = shapely.geometry.GeometryCollection([multipolygon1,line2,pnt1])

geo_collection

Shapely μέσω shape¶

Μπορούμε να δημιουργήσουμε αντικείμενα shapely μέσω της συνάρτησης shapely.geometry.shape η οποία δέχεται σαν όρισμα ένα λεξικό μορφής GEOJSON το οποίο πρέπει να έχει τις εξής δύο ιδιότητες.

Την ιδιότητα "type" που περιγράφει τον τύπο γεωμετρίας
Την ιδιότητα "coordinates" που περιγράφει τις γεωμετρίες και τις συντεταγμένες τους σαν λίστες ή πλειάδες.

d = {"type": "Point", "coordinates": (0, 1)}
shapely.geometry.shape(d)

Δημιουργία αντικειμένου MultiPolygon μέσω GEOJSON λεξικού:

d = {
  "type": "MultiPolygon",
  "coordinates": [
    [
      [[40, 40], [20, 45], [45, 30], [40, 40]]
    ],
    [
      [[20, 35], [10, 30], [10, 10], [30, 5], [45, 20], 
      [20, 35]], 
      [[30, 20], [20, 15], [20, 25], [30, 20]]
    ]
  ]
}
pol3 = shapely.geometry.shape(d)
pol3

Γεωμετρικός τύπος δεδομένων¶

H ιδιότητα .geom_type ενός αντικειμένου shapely περιγράφει τον γεωμετρικό τύπο της:

pol1.geom_type

'Polygon'

line1.geom_type

'LineString'

geo_collection.geom_type

'GeometryCollection'

Συντεταγμένες¶

Για να ανακτήσουμε τις συντεταγμένες ενός shapely αντικειμένου καλούμε με διαφορετικό τρόπο τις ανάλογες συναρτήσεις ανάλογα την πολυπλοκότητα του κάθε τύπου.

Για ένα Point object καλούμε άμεσα την ιδιότητα coords η οποία επιστρέφει ένα shapely.coords.CoordinateSequence object.

pnt1.coords

<shapely.coords.CoordinateSequence at 0x7f970f4f2ad0>

Μπορούμε να λάβουμε ως λίστα τις πλειάδες συντεταγμένων που το απαρτίζουν

list(pnt1.coords)

[(2.0, 0.5)]

Αντίστοιχα για γραμμή

list(line1.coords)

[(2.0, 0.5), (1.0, 1.0), (-1.0, 0.0), (1.0, 0.0)]

Για να ελέγξουμε σε ένα αντικείμενο συλλογών γεωμετρίας (MultiPoint, MultiPolygon κτλ) πόσες επιμέρους γεωμετρίες περιλαμβάνει:

len(line2.geoms)

line2

Για να λάβουμε το πρώτο αντικείμενο γεωμετρίας:

line2.geoms[0]

type(line2.geoms[0])

shapely.geometry.linestring.LineString

line2.geoms[0].geom_type

'LineString'

Και λαμβάνουμε τις συντεταγμένες της πρώτης γεωμετρίας:

list(line2.geoms[0].coords)

[(2.0, 0.5), (1.0, 1.0), (-1.0, 0.0), (1.0, -1.0)]

ή της δεύτερης

list(line2.geoms[1].coords)

[(-2.0, 1.0), (2.0, -0.5), (3.0, 1.0)]

Για τα πολύγωνα ακολουθούμε διαφορετική προσέγγιση. Ένα πολύγωνο αποτελεί από το εξωτερικό περίγραμμα (exterior) ή και ένα ή περισσότερα περιγράμματα εσωτερικών τρυπών (interiors). Κατά συνέπεια έχουμε συντεταγμένες που περιγράφουν το κάθε περίγραμμα.

Το εξωτερικό περίγραμμα του pol1 αντικειμένου:

pol1.exterior

Και οι συντεταγμένες του:

pol1

list(pol1.exterior.coords)

[(0.0, 0.0), (0.0, -1.0), (7.5, -1.0), (7.5, 0.0), (0.0, 0.0)]

Το pol1 όμως δεν έχει εσωτερικές τρύπες γι αυτό και το παρακάτω επιστρέφει μηδέν.

len(pol1.interiors)

Έστω το παρακάτω MultiPolygon object που δημιουργήσαμε σε προηγούμενο στάδιο.

pol3

Περιλαμβάνει δύο ξεχωριστά γεωμετρικά αντικείμενα:

len(pol3.geoms)

Ας πάρουμε το πρώτο:

pol3.geoms[0]

Δεν περιλαμβάνει καμία τρύπα στο εσωτερικό του. Ας πάρουμε τις συντεταγμένες από το περίγραμμά του (exterior)

pol3.geoms[0].exterior.coords

<shapely.coords.CoordinateSequence at 0x7f9750d313c0>

Και ας τις επιστρέψουμε σαν λίστα πλειάδων από ζεύγη της μορφής (x,y)

list(pol3.geoms[0].exterior.coords)

[(40.0, 40.0), (20.0, 45.0), (45.0, 30.0), (40.0, 40.0)]

Ας δοκιμάσουμε το δεύτερο αντικείμενο γεωμετρίας. Ας το δούμε:

pol3.geoms[1]

Λαμβάνουμε τις συντεταγμένες του εξωτερικού περιγράμματος:

list(pol3.geoms[1].exterior.coords)

[(20.0, 35.0),
 (10.0, 30.0),
 (10.0, 10.0),
 (30.0, 5.0),
 (45.0, 20.0),
 (20.0, 35.0)]

Ας δούμε πόσες τρύπες έχει στο εσωτερικό του:

len(pol3.geoms[1].interiors)

Παίρνουμε το περίγραμμα της τρύπας

pol3.geoms[1].interiors[0]

Και τις συντεταγμένες της

list(pol3.geoms[1].interiors[0].coords)

[(30.0, 20.0), (20.0, 15.0), (20.0, 25.0), (30.0, 20.0)]

Ιδιότητες αντικειμένων¶

Υπολογισμός ορίων (bounds)¶

geo_collection

geo_collection.bounds

(-2.0, -1.0, 45.0, 45.0)

shapely.geometry.box(*geo_collection.bounds)

line1

line1.bounds

(-1.0, 0.0, 2.0, 1.0)

pnt1

pnt1.bounds

(2.0, 0.5, 2.0, 0.5)

list(pnt1.coords)

[(2.0, 0.5)]

Υπολογισμός μήκους γραμμής¶

line1

line1.length

5.354101966249685

line2.geoms[0]

line2.geoms[0].length

5.5901699437494745

Υπολογισμός εμβαδού¶

pol1.area

7.5

pol2.area

6.25

Νέες γεωμετρίες¶

Κεντροειδές πολυγώνου (centroid)¶

pol2

pol2.centroid

shapely.geometry.GeometryCollection([pol2, pol2.centroid])

Περιμετρική ζώνη (buffer)¶

pnt1.buffer(5)

shapely.geometry.GeometryCollection([pnt1,pnt1.buffer(5)])

pol1.buffer(5)

shapely.geometry.GeometryCollection([pol1,pol1.buffer(5)])

Convex hull¶

pol3

pol3.convex_hull

Σχέσεις μεταξύ αντικειμένων¶

shapely.geometry.GeometryCollection([pol1,pol3])

pol3

pol1.intersects(pol3)

False

pol1.intersects(pol2)

True

shapely.geometry.GeometryCollection([pol1,pol2])

Γεωμετρικές πράξεις¶

x = shapely.geometry.Point((0, 0)).buffer(1)
y = shapely.geometry.Point((1, 0)).buffer(1)
shapely.geometry.GeometryCollection([x, y])

x.intersection(y)

x.difference(y)

x.union(y)

pol1.union(pol2)

Υπολογισμός απόστασης ανάμεσα σε δύο αντικείμενα

shapely.geometry.GeometryCollection([pol1, pol3])

pol1.distance(pol3)

10.307764064044152

pol3.distance(pol1) 

10.307764064044152

import pyproj

from shapely.geometry import Point
from shapely.ops import transform

wgs84_pt = Point(39.35858398397631, 22.932070043920394)

wgs84 = pyproj.CRS('EPSG:4326')
greek_grid = pyproj.CRS('EPSG:2100')

project = pyproj.Transformer.from_crs(wgs84, greek_grid, always_xy=True).transform
projected_point = transform(project, wgs84_pt)

print(projected_point)

POINT (2087909.8290560723 2620022.621513317)

list(wgs84_pt.coords)[0]

(39.35858398397631, 22.932070043920394)

import folium

coords = list(wgs84_pt.coords)[0]

#Create the map
my_map = folium.Map(location = coords, zoom_start = 13)

# Add marker
folium.Marker(coords, popup = 'Volos').add_to(my_map)

#Display the map
my_map

Make this Notebook Trusted to load map: File -> Trust Notebook

Η βιβλιοθήκη Geopandas¶

Μέχρι τώρα είδαμε πως μπορούμε να διαχειριζόμαστε γεωμετρικά δεδομένα με την βιβλιοθήκη shapely.

Όμως τα γεωγραφικά δεδομένα δεν περιλαμβάνουν μόνο της γεωγραφική πληροφορία και την γεωμετρική τους δομή αλλά συνοδεύονται από μια σειρά περιγραφικών δεδομένων.

Η βιβλιοθήκη geopandas είναι μια βιβλιοθήκη της Python η οποία υποστηρίζει την ανάγνωση, επεξεργασία, ανάλυση και εγγραφή γεωγραφικών και παράλληλα περιγραφικών δεδομένων. Αποτελεί επέκταση της βιβλιοθήκης pandas και «κληρονομεί» τα χαρακτηριστικά και τις δυνατότητές της.

Στην υφιστάμενη δομή της pandas, η geopandas υποστηρίζει την γεωμετρία με την προσθήκη μιας νέας στήλης και το γεωγραφικό σύστημα αναφοράς.

alt text

Εισάγουμε την βιβλιοθήκη με τον παρακάτω τρόπο

import geopandas as gpd

Για να διαβάσουμε ένα αρχείο shapefile χρησιμοποιούμε την συνάρτηση gpd.read_file. Το αντικείμενο που προκύπτει είναι τύπου geopandas.geodataframe.GeoDataFrame

# Import shapefile using geopandas
dhmoi = gpd.read_file("../docs/dhmoi.gpkg")

type(dhmoi)

geopandas.geodataframe.GeoDataFrame

Η στήλη της γεωμετρίας ονομάζεται προκαθορισμένα geometry και είναι αντικείμενο geopandas.geoseries.GeoSeries που περιλαμβάνει αντικείμενα shapely.

type(dhmoi["geometry"])

geopandas.geoseries.GeoSeries

Με την μέθοδο geom_type μπορούμε να δούμε τον γεωμετρικό τύπο κάθε εγγραφής.

dhmoi.geom_type

    MultiPolygon
    MultiPolygon
    MultiPolygon
    MultiPolygon
    MultiPolygon
           ...     
  MultiPolygon
  MultiPolygon
  MultiPolygon
  MultiPolygon
  MultiPolygon
Length: 325, dtype: object

Ας δούμε τις αρχικές εγγραφές του αρχείου

dhmoi.head()

	OBJECTID	X	Y	Name	CodeELSTAT	PopM01	PopF01	PopTot01	UnemrM01	UnemrF01	UnemrT01	PrSect01	Foreig01	Income01	Perif	geometry
0	1	616259.007833	4.551127e+06	KOMOTINIS	0101	29967.0	31534.0	61501.0	4.726948	6.051873	5.221281	21.817852	1.996780	10969.315256	1	MULTIPOLYGON (((631058.735 4574912.809, 631431...
1	2	644783.135624	4.561364e+06	ARRIANON	0102	8952.0	9307.0	18259.0	1.955813	2.390181	2.154240	92.912436	0.038337	6398.641486	1	MULTIPOLYGON (((661976.597 4554253.076, 662118...
2	3	602100.950461	4.555689e+06	IASMOU	0103	7290.0	7561.0	14851.0	5.653884	6.671554	6.062390	70.206767	1.339977	6608.897469	1	MULTIPOLYGON (((603623.489 4566376.502, 603657...
3	4	633712.468442	4.539548e+06	MARONEIAS - SAPON	0104	8284.0	8342.0	16626.0	5.840957	10.270499	7.421934	57.748085	0.890172	6700.375345	1	MULTIPOLYGON (((650998.525 4557469.860, 651446...
4	5	517831.195188	4.572544e+06	DRAMAS	0201	28041.0	29326.0	57367.0	8.831689	14.176926	10.851115	7.298838	2.961633	11299.106498	2	MULTIPOLYGON (((521986.700 4601855.354, 522271...

dhmoi.plot()

<AxesSubplot:>

Χαρτογραφική απόδοση σε διαδραστικό χάρτη

dhmoi.explore(legend=False)

Make this Notebook Trusted to load map: File -> Trust Notebook

Εκτύπωση λεπτομερειών για το γεωγραφικό σύστημα αναφοράς το οποίο όπως φαίνεται παρακάτω είναι το ΕΓΣΑ “87 (EPSG:2100)

dhmoi.crs

<Derived Projected CRS: PROJCS["unknown",GEOGCS["unknown",DATUM["Greek_Geo ...>
Name: unknown
Axis Info [cartesian]:
- [east]: Easting (metre)
- [north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unnamed
- method: Transverse Mercator
Datum: Greek Geodetic Reference System 1987
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

Μπορούμε να μετασχηματίσουμε τα δεδομένα σε ένα διαφορετικό προβολικό σύστημα με κλήση της μεθόδου to_crs

dhmoi_wgs84 = dhmoi.to_crs(4326)
dhmoi_wgs84.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

dhmoi.geom_type

    MultiPolygon
    MultiPolygon
    MultiPolygon
    MultiPolygon
    MultiPolygon
           ...     
  MultiPolygon
  MultiPolygon
  MultiPolygon
  MultiPolygon
  MultiPolygon
Length: 325, dtype: object

Εκτύπωση των γεωγραφικών ορίων του αρχείου oikismoi που είναι στο σύστημα αναφορά ΕΓΣΑ “87.

dhmoi.total_bounds

array([ 104040.7266, 3850938.    , 1007943.    , 4624010.2508])

Εκτύπωση των γεωγραφικών ορίων του αρχείου oikismoi που είναι στο σύστημα αναφορά WGS’84.

dhmoi_wgs84.total_bounds

array([19.37377196, 34.80060523, 29.64123782, 41.74911332])

Η μέθοδος shape μας επιστρέφει τον πλήθος των γραμμών (325) και των στηλών (15).

dhmoi.shape

(325, 16)

Μπορούμε να εγγράψουμε ένα Geopandas Dataframe

dhmoi.to_file("dhmoi.geojson", driver="GeoJSON")

/home/leonidas/anaconda3/envs/book/lib/python3.10/site-packages/geopandas/io/file.py:362: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.
  pd.Int64Index,

dhmoi.sort_values(by='PopTot01', ascending=False)

	OBJECTID	X	Y	Name	CodeELSTAT	PopM01	PopF01	PopTot01	UnemrM01	UnemrF01	UnemrT01	PrSect01	Foreig01	Income01	Perif	geometry
192	193	476676.870086	4.204338e+06	ATHINAION	4501	374900.0	414266.0	789166.0	4.890244	6.640842	5.624998	0.515649	17.426785	14497.393539	45	MULTIPOLYGON (((478164.675 4208572.203, 478015...
22	23	411331.939920	4.497210e+06	THESSALONIKIS	0701	185771.0	211385.0	397156.0	5.402357	7.615427	6.336333	0.834415	7.200951	13388.724605	7	MULTIPOLYGON (((408112.992 4500238.217, 408264...
147	148	306893.782874	4.231792e+06	PATREON	3701	104429.0	106065.0	210494.0	7.613916	9.985632	8.493895	2.567481	5.301864	12699.732926	37	MULTIPOLYGON (((312258.406 4245213.000, 312382...
245	246	468574.034133	4.199691e+06	PEIRAIOS	5101	87362.0	94571.0	181933.0	5.742561	7.683770	6.480184	0.491722	8.515141	13526.486840	51	MULTIPOLYGON (((469779.406 4198249.000, 469760...
301	302	599279.522964	3.900401e+06	IRAKLEIOU	7101	81067.0	82048.0	163115.0	5.047587	8.832661	6.557167	7.477305	4.048727	12479.911946	71	MULTIPOLYGON (((602945.188 3911414.000, 602942...
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
299	300	824000.991323	4.016217e+06	CHALKIS	6905	165.0	130.0	295.0	1.265823	4.166667	1.941748	31.683168	17.966102	8888.138597	69	MULTIPOLYGON (((830115.688 4011995.000, 830198...
270	271	659987.832319	4.025277e+06	ANAFIS	6002	135.0	137.0	272.0	1.369863	2.941176	1.869159	43.809524	1.102941	9062.404166	60	MULTIPOLYGON (((669453.375 4014893.000, 669444...
272	273	599675.697252	4.059438e+06	SIKINOU	6004	118.0	120.0	238.0	0.000000	4.545455	1.333333	27.027027	1.680672	7566.666874	60	MULTIPOLYGON (((593612.000 4054459.000, 593596...
275	276	762387.857032	4.149989e+06	AGATHONISIOU	6102	92.0	60.0	152.0	12.765957	10.000000	12.280702	56.000000	9.210526	5421.733909	61	MULTIPOLYGON (((763781.625 4146430.000, 763859...
320	321	507562.964610	3.855768e+06	GAVDOU	7403	49.0	32.0	81.0	0.000000	0.000000	0.350000	51.351351	9.876543	8498.080328	74	MULTIPOLYGON (((505925.688 3858625.000, 506062...

325 rows × 16 columns

Ορισμός του index στην στήλη «CodeELSTAT»

dhmoi = dhmoi.set_index("CodeELSTAT")

dhmoi

	OBJECTID	X	Y	Name	PopM01	PopF01	PopTot01	UnemrM01	UnemrF01	UnemrT01	PrSect01	Foreig01	Income01	Perif	geometry
CodeELSTAT
0101	1	616259.007833	4.551127e+06	KOMOTINIS	29967.0	31534.0	61501.0	4.726948	6.051873	5.221281	21.817852	1.996780	10969.315256	1	MULTIPOLYGON (((631058.735 4574912.809, 631431...
0102	2	644783.135624	4.561364e+06	ARRIANON	8952.0	9307.0	18259.0	1.955813	2.390181	2.154240	92.912436	0.038337	6398.641486	1	MULTIPOLYGON (((661976.597 4554253.076, 662118...
0103	3	602100.950461	4.555689e+06	IASMOU	7290.0	7561.0	14851.0	5.653884	6.671554	6.062390	70.206767	1.339977	6608.897469	1	MULTIPOLYGON (((603623.489 4566376.502, 603657...
0104	4	633712.468442	4.539548e+06	MARONEIAS - SAPON	8284.0	8342.0	16626.0	5.840957	10.270499	7.421934	57.748085	0.890172	6700.375345	1	MULTIPOLYGON (((650998.525 4557469.860, 651446...
0201	5	517831.195188	4.572544e+06	DRAMAS	28041.0	29326.0	57367.0	8.831689	14.176926	10.851115	7.298838	2.961633	11299.106498	2	MULTIPOLYGON (((521986.700 4601855.354, 522271...
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
7403	321	507562.964610	3.855768e+06	GAVDOU	49.0	32.0	81.0	0.000000	0.000000	0.350000	51.351351	9.876543	8498.080328	74	MULTIPOLYGON (((505925.688 3858625.000, 506062...
7404	322	475420.582085	3.905918e+06	KANTANOU - SELINOU	3271.0	3031.0	6302.0	2.312500	2.234637	2.284569	61.525841	14.680209	7989.958401	74	MULTIPOLYGON (((469969.812 3897335.000, 469955...
7405	323	465887.990746	3.920087e+06	KISSAMOU	5925.0	5545.0	11470.0	4.847458	5.894591	5.191257	45.701249	7.916303	8750.030129	74	MULTIPOLYGON (((457703.500 3902963.000, 457852...
7406	324	484468.969583	3.924051e+06	PLATANIA	9217.0	8647.0	17864.0	4.437401	8.586262	5.938448	51.305684	8.889386	8120.027122	74	MULTIPOLYGON (((476162.812 3949684.000, 476253...
7407	325	508083.904768	3.901847e+06	SFAKION	1288.0	1131.0	2419.0	1.920439	2.839117	2.198853	63.929619	4.960728	6602.391386	74	MULTIPOLYGON (((502587.434 3912125.553, 502774...

325 rows × 15 columns

Υπολογισμός έκτασης (σε \(km^2\)) σε μία νέα στήλη με το όνομα area

dhmoi["area"] = dhmoi.area*10e-6

dhmoi

	OBJECTID	X	Y	Name	PopM01	PopF01	PopTot01	UnemrM01	UnemrF01	UnemrT01	PrSect01	Foreig01	Income01	Perif	geometry	area
CodeELSTAT
0101	1	616259.007833	4.551127e+06	KOMOTINIS	29967.0	31534.0	61501.0	4.726948	6.051873	5.221281	21.817852	1.996780	10969.315256	1	MULTIPOLYGON (((631058.735 4574912.809, 631431...	6470.182562
0102	2	644783.135624	4.561364e+06	ARRIANON	8952.0	9307.0	18259.0	1.955813	2.390181	2.154240	92.912436	0.038337	6398.641486	1	MULTIPOLYGON (((661976.597 4554253.076, 662118...	7724.393191
0103	3	602100.950461	4.555689e+06	IASMOU	7290.0	7561.0	14851.0	5.653884	6.671554	6.062390	70.206767	1.339977	6608.897469	1	MULTIPOLYGON (((603623.489 4566376.502, 603657...	4870.493148
0104	4	633712.468442	4.539548e+06	MARONEIAS - SAPON	8284.0	8342.0	16626.0	5.840957	10.270499	7.421934	57.748085	0.890172	6700.375345	1	MULTIPOLYGON (((650998.525 4557469.860, 651446...	6441.606153
0201	5	517831.195188	4.572544e+06	DRAMAS	28041.0	29326.0	57367.0	8.831689	14.176926	10.851115	7.298838	2.961633	11299.106498	2	MULTIPOLYGON (((521986.700 4601855.354, 522271...	8401.321201
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
7403	321	507562.964610	3.855768e+06	GAVDOU	49.0	32.0	81.0	0.000000	0.000000	0.350000	51.351351	9.876543	8498.080328	74	MULTIPOLYGON (((505925.688 3858625.000, 506062...	343.773091
7404	322	475420.582085	3.905918e+06	KANTANOU - SELINOU	3271.0	3031.0	6302.0	2.312500	2.234637	2.284569	61.525841	14.680209	7989.958401	74	MULTIPOLYGON (((469969.812 3897335.000, 469955...	3746.368061
7405	323	465887.990746	3.920087e+06	KISSAMOU	5925.0	5545.0	11470.0	4.847458	5.894591	5.191257	45.701249	7.916303	8750.030129	74	MULTIPOLYGON (((457703.500 3902963.000, 457852...	3421.398731
7406	324	484468.969583	3.924051e+06	PLATANIA	9217.0	8647.0	17864.0	4.437401	8.586262	5.938448	51.305684	8.889386	8120.027122	74	MULTIPOLYGON (((476162.812 3949684.000, 476253...	4925.063676
7407	325	508083.904768	3.901847e+06	SFAKION	1288.0	1131.0	2419.0	1.920439	2.839117	2.198853	63.929619	4.960728	6602.391386	74	MULTIPOLYGON (((502587.434 3912125.553, 502774...	4667.779146

325 rows × 16 columns

Υπολογισμός του centroid κάθε δήμου σε μια νέα στήλη (centroid)

dhmoi['centroid'] = dhmoi.centroid

dhmoi

	OBJECTID	X	Y	Name	PopM01	PopF01	PopTot01	UnemrM01	UnemrF01	UnemrT01	PrSect01	Foreig01	Income01	Perif	geometry	area	centroid
CodeELSTAT
0101	1	616259.007833	4.551127e+06	KOMOTINIS	29967.0	31534.0	61501.0	4.726948	6.051873	5.221281	21.817852	1.996780	10969.315256	1	MULTIPOLYGON (((631058.735 4574912.809, 631431...	6470.182562	POINT (616259.008 4551127.241)
0102	2	644783.135624	4.561364e+06	ARRIANON	8952.0	9307.0	18259.0	1.955813	2.390181	2.154240	92.912436	0.038337	6398.641486	1	MULTIPOLYGON (((661976.597 4554253.076, 662118...	7724.393191	POINT (644783.136 4561364.036)
0103	3	602100.950461	4.555689e+06	IASMOU	7290.0	7561.0	14851.0	5.653884	6.671554	6.062390	70.206767	1.339977	6608.897469	1	MULTIPOLYGON (((603623.489 4566376.502, 603657...	4870.493148	POINT (602100.950 4555689.138)
0104	4	633712.468442	4.539548e+06	MARONEIAS - SAPON	8284.0	8342.0	16626.0	5.840957	10.270499	7.421934	57.748085	0.890172	6700.375345	1	MULTIPOLYGON (((650998.525 4557469.860, 651446...	6441.606153	POINT (633712.468 4539548.490)
0201	5	517831.195188	4.572544e+06	DRAMAS	28041.0	29326.0	57367.0	8.831689	14.176926	10.851115	7.298838	2.961633	11299.106498	2	MULTIPOLYGON (((521986.700 4601855.354, 522271...	8401.321201	POINT (517831.195 4572544.478)
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
7403	321	507562.964610	3.855768e+06	GAVDOU	49.0	32.0	81.0	0.000000	0.000000	0.350000	51.351351	9.876543	8498.080328	74	MULTIPOLYGON (((505925.688 3858625.000, 506062...	343.773091	POINT (507562.965 3855767.657)
7404	322	475420.582085	3.905918e+06	KANTANOU - SELINOU	3271.0	3031.0	6302.0	2.312500	2.234637	2.284569	61.525841	14.680209	7989.958401	74	MULTIPOLYGON (((469969.812 3897335.000, 469955...	3746.368061	POINT (475420.582 3905917.949)
7405	323	465887.990746	3.920087e+06	KISSAMOU	5925.0	5545.0	11470.0	4.847458	5.894591	5.191257	45.701249	7.916303	8750.030129	74	MULTIPOLYGON (((457703.500 3902963.000, 457852...	3421.398731	POINT (465887.991 3920086.947)
7406	324	484468.969583	3.924051e+06	PLATANIA	9217.0	8647.0	17864.0	4.437401	8.586262	5.938448	51.305684	8.889386	8120.027122	74	MULTIPOLYGON (((476162.812 3949684.000, 476253...	4925.063676	POINT (484468.970 3924050.753)
7407	325	508083.904768	3.901847e+06	SFAKION	1288.0	1131.0	2419.0	1.920439	2.839117	2.198853	63.929619	4.960728	6602.391386	74	MULTIPOLYGON (((502587.434 3912125.553, 502774...	4667.779146	POINT (508083.905 3901847.180)

325 rows × 17 columns

Χαρτογραφική απόδοση του δείκτη ανεργίας (στήλη UnemrT01) ανά δήμο

dhmoi.plot("UnemrT01", legend=True)

<AxesSubplot:>

 dhmoi.plot("Income01", scheme='quantiles', cmap='YlOrRd',  legend=True, figsize=(12, 10))

<AxesSubplot:>

Μπορούμε να οπτικοποιήσουμε διαφορετική στήλη τύπου geopandas.geoseries.GeoSeries αφού πρώτα την ορίσουμε με την μέθοδο set_geometry.

dhmoi = dhmoi.set_geometry("centroid")
dhmoi.plot("UnemrT01", legend=True)

<AxesSubplot:>

Ορισμός περιμετρικής ζώνης διαμέτρου 20χλμ γύρω από κάθε centroid.

dhmoi = dhmoi.set_geometry("centroid") # ορισμός default γεωμετρίας η στήλη centroid
dhmoi["buffered"] = dhmoi.buffer(20000) #εκτέλεση περιμετρικών ζωνών

dhmoi = dhmoi.set_geometry("buffered") # ορισμός default γεωμετρίας η στήλη buffered
dhmoi.plot(legend=True) # οπτικοποίηση

<AxesSubplot:>

Ξανα ορίζουμε σαν προκαθορισμένη στήλη γεωμετρία την στήλη geometry.

dhmoi = dhmoi.set_geometry("geometry")

Μπορούμε να φιλτράρουμε γραμμές και στήλες όπως και στην pandas χρησιμοποιώντας το index και το όνομα στήλης:

lesvos = dhmoi.loc["5301", "geometry"]

lesvos

type(lesvos)

shapely.geometry.multipolygon.MultiPolygon

Η παραπάνω διαδικασία επιλογής επέστρεψε ένα shapely MultiPolygon object όπως φαίνεται.

Με την μέθοδο dissolve μπορούμε να συγχωνεύσουμε γεωμετρικά αντικείμενα βάση κάποιας στήλης. Αντικείμενα με ίδια τιμή θα συγχωνευθούν σε ενιαίο γεωμετρικό αντικείμενο.

dhmoi.head()

	OBJECTID	X	Y	Name	PopM01	PopF01	PopTot01	UnemrM01	UnemrF01	UnemrT01	PrSect01	Foreig01	Income01	Perif	geometry	area	centroid	buffered
CodeELSTAT
0101	1	616259.007833	4.551127e+06	KOMOTINIS	29967.0	31534.0	61501.0	4.726948	6.051873	5.221281	21.817852	1.996780	10969.315256	1	MULTIPOLYGON (((631058.735 4574912.809, 631431...	6470.182562	POINT (616259.008 4551127.241)	POLYGON ((636259.008 4551127.241, 636162.702 4...
0102	2	644783.135624	4.561364e+06	ARRIANON	8952.0	9307.0	18259.0	1.955813	2.390181	2.154240	92.912436	0.038337	6398.641486	1	MULTIPOLYGON (((661976.597 4554253.076, 662118...	7724.393191	POINT (644783.136 4561364.036)	POLYGON ((664783.136 4561364.036, 664686.830 4...
0103	3	602100.950461	4.555689e+06	IASMOU	7290.0	7561.0	14851.0	5.653884	6.671554	6.062390	70.206767	1.339977	6608.897469	1	MULTIPOLYGON (((603623.489 4566376.502, 603657...	4870.493148	POINT (602100.950 4555689.138)	POLYGON ((622100.950 4555689.138, 622004.645 4...
0104	4	633712.468442	4.539548e+06	MARONEIAS - SAPON	8284.0	8342.0	16626.0	5.840957	10.270499	7.421934	57.748085	0.890172	6700.375345	1	MULTIPOLYGON (((650998.525 4557469.860, 651446...	6441.606153	POINT (633712.468 4539548.490)	POLYGON ((653712.468 4539548.490, 653616.163 4...
0201	5	517831.195188	4.572544e+06	DRAMAS	28041.0	29326.0	57367.0	8.831689	14.176926	10.851115	7.298838	2.961633	11299.106498	2	MULTIPOLYGON (((521986.700 4601855.354, 522271...	8401.321201	POINT (517831.195 4572544.478)	POLYGON ((537831.195 4572544.478, 537734.890 4...

nomoi = dhmoi.dissolve(by='Perif')

dhmoi.plot()

<AxesSubplot:>

nomoi.plot()

<AxesSubplot:>

Επιπλεόν κατά την συγχώνευση μπορούμε να εφαρμόσουμε μια συνάρτηση στις στήλες π.χ. mean ή sum.

nomoi = dhmoi.dissolve(by='Perif', aggfunc='mean')

nomoi.plot(column = 'Income01', scheme='quantiles', cmap='YlOrRd');

Βιβλιογραφία:¶

Kalogirou S., 2020, Local Correlation, Spatial Inequalities, Geographically Weighted Regression and Other Tools (R package)
GeoPandas, https://geopandas.org/, Πρόσβαση: 29/05/2022
Prapas I., Analyze Geospatial Data in Python: GeoPandas and Shapely, https://www.learndatasci.com/tutorials/geospatial-data-python-geopandas-shapely/, Πρόσβαση: 29/05/2022
Dorman M., 2022, VECTOR LAYERS, Geometries (shapely), https://geobgu.xyz/py/shapely.html, Πρόσβαση: 29/05/2022
Wasser, Leah, Korinek, Nathan, & Palomino, Jenny. (2021). earthlab/earth-analytics-intermediate-earth-data-science-textbook: one more license fix (1.0.4) [Computer software]. Zenodo. https://doi.org/10.5281/zenodo.5571001

Προγραμματισμός Η/Υ

9. Γεωεπεξεργασία διανυσματικών δεδομένων

Περιεχόμενα

9. Γεωεπεξεργασία διανυσματικών δεδομένων¶

Python βιβλιοθήκες για διανυσματικά δεδομένα¶

Η βιβλιοθήκη Shapely¶

shapely from WKT¶

Πολύγωνα¶

Σημεία¶

Γραμμές¶

Shapely μέσω συναρτήσεων¶

Shapely μέσω shape¶

Γεωμετρικός τύπος δεδομένων¶

Συντεταγμένες¶

Ιδιότητες αντικειμένων¶

Υπολογισμός ορίων (bounds)¶

Υπολογισμός μήκους γραμμής¶

Υπολογισμός εμβαδού¶

Νέες γεωμετρίες¶

Κεντροειδές πολυγώνου (centroid)¶

Περιμετρική ζώνη (buffer)¶

Convex hull¶

Σχέσεις μεταξύ αντικειμένων¶

Γεωμετρικές πράξεις¶

Η βιβλιοθήκη Geopandas¶

Βιβλιογραφία:¶